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FAST FOURIER TRANSFORM (FFT) 


DECUS Program Library Write-up DECUS No. 8-250 

OPERATING INSTRUCTIONS 
Starting Address = 0200 

Program immediately halts at location 0200. Set the switch register for desired mode of 
operation. 

Switch settings 

SWO up-Convert and sample on RC timing clock flag. 

SWO down- Convert and sample every 0.01 seconds. 

N many samples are taken and then displayed on the scope. If 

SWll up—Keep displaying input time series until lowered. 

SW11 down- Display only once. 

Compute the frequency spectrum of the input, in less than 8 seconds, and display the output 
on the scope. 

SW10 up-Keep displaying until lowered. 

SW10 down-Display once then return to beginning of the program. 

SW9 up-Display histogram. 

SW9 down--Display only points. 

Knob 34-Controls 8 step vertical scaling. 

Knob 35—-Controls 32 step horizontal scaling. 

Knob 36-Controls 1024 step starting frequency value. 

The Fast Fourier Transform 

Very often it is necessary to know more about a signal than just what it or n-many averaged out 
look like. For purposes of analysis, it is often helpful to know what the frequency components 
of the waveform are. A good example of this is the analysis of electroencephalograms (EEG's) to 
determine sleep stages. Determining the frequency components of a wave involves a shift from 
the time domain into the frequency domain. This can be done using Fourier analysis. 

The underlying assumption of the whole field of Fourier analysis is that every piecewise continuous, 
absolutely integrable function can be approximated by an infinite series of sine and cosine waves, or; 
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n=0 


(equ. 1 .) 
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In which T is the period of the function, and Wq, and Wq-2/1 . 

By noting the relationship between the trigonometric functions and the complex exponentials; 


. 1 * = 


cos(x) + jsin(x) 


sin(wt) = (e’ Wt — e~' Wt ) 

cos(wt) = y (e' wt + e“ iwf ) (equ. 2 .) 

where; j= 

and e =2.72 is the base for natural logarithms, 
the Fourier series (equ. 1. ) becomes; 
n=0O 

f(0 = n e inw 0 f 


(equ. 3.) 


where; 


n=0 


c = 1/2 (a - b ) = 
n n n 


jJ0 

e^n 


- l/2(o n 2 + b 2 ) '/ 2 

n n 


0 = arctan(-b /a ). 
n n n 

If we go to the limiting cases, and let T approach infinity, the Fourier series becomes the 
Fourier integral. This integration process is known as the Fourier Transform and is defined as; 

00 

fW - F(f(t) ) =Cf(t) e'' wt dt (equ. 4.). 
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Equation 4, as was mentioned above, was for piecewise continuous functions. However, if 
digital analysis techniques are to be used in the analysis of the waveform, then it is necessary 
that the data be sampled in order to produce a time series of discrete samples which can be fed 
into a digital computer. The fundamental theorem of sampling theory says that if a waveform 
is frequency band limited and the samples are taken at a rate at least twice that of the highest 
frequency present, then the time series completely and accurately represents the continuous 
waveform. These samples need not be equally spaced in time, but for the purposes of analysis, 
and simplicity, they will be considered as if they were equally spaced. The frequency twice 
the highest present in the waveform is known as the Nyquist or aliasing frequency. The Fourier 
transform that is taken of this time series is known as the Discrete Fourier Transform (DFT). 


The formation of the time series may be looked upon as the multiplication of a continuous 
function, x(t), by a sampling function, s(t). In this case, s(t) is actually an infinite series of 
impulses. Mathematically, if two functions are multiplied, then their Fourier transforms are 
subjected to a process called convolution. This makes the frequency function, F(x(t) ) = x(w), 

• a periodic function, where each period contains complete information about the frequency 

components of the original continuous waveform. This frequency spectrum has a period equal to 
the inverse of the spacing between samples. Because of this periodicity, only one period need 
be calculated. Graphical representation of this is provided in figure 1. 

The frequency components of the DFT are defined by equ. 5. 


k=N-l 

S(n) =At ^ X(k) e = i 2 ‘I r ( nk )/N 
k=0 


n=0,1,2,..., N-l 
(equ. 5.) 


Where X(k) is the K time sample, and N is the total number of samples taken. 


If T is taken to be the time interval between the samples, then the true frequency is given by 
equ. 6. 7 

1 •“ k /<NT) 6.) 

Thus it is very important to know both N and T if the interpretation of S(n) is to be meaningful. 
Equation 5. is often written in the form shown by; 

k= N-l 

s(n) = a * I — x, w rk 

k=0 k 

where; W = e"i 21 T /N . 
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Evaluation of equ. 7. would seem to imply N complex multiplications. However, Cooley and 
Tukey devised a method of evaluating it if N is an integer power of 2. This process involves no 
more than Nlog 2 N complex multiplications. Because of the vast reduction in computation time, 
this algorithm became known as the Fast Fourier Transform (FFT). Computationally, equ. 7. is # 
nothing more than a complex matrix multiplication of the form; 


r=0,1,2,.,., N-l 

(equ. 7) 
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(equ. 8.) 


H 

where; 


and; 


■ M M 

[ s< n>] 


M is 


^X Q (k)J are N 


an N*N square matrix. 


column matrices. 


The term^t is only a scalar multiplier and is of little computational value, and thus, ignored. 
It should also be noted that the elements of [S(n)] will be in binary bit reverse order if 
IX M] is properly ordered, and vice versa. We see that since all the elements of the 
matrfces are complex, a matrix multiply routine would take up too much space and too much 
time . So the implementation of the FFT will be discussed. 


Start by drawing a vertical array of (=2^) nodes. This array is called X , so that sample X, 
is associated with the K node of X . The number k will be expressed in binary form, and will 
take on values ranging from 0 to N-1. Now construct g many of these arrays to the right of 
the original. Each node has a number within it, and a j^umber associated with it (but not inside 
it). The number X, is the number associated with the k node in the e array, (e is now just 
an index, and not to be confused with the natural logarithm base.) 


Looking at figure 2., a tree-graph, every node has a solid and a dashed line drawn to it from two 
different nodes in the previous array. The solid lines stand for multiplication, and the dashed 
ones sta£jd for addition,^jTie two processes that are involved in a complex multiply. If the number 
in the k node of the e array is C, and the solid line is drawn to that node from node A, and 
the dashed line comes from node B, then X^ is formed by taking the number associated with A 
and multiplying it by C, and then added to the number associated with B. Letting lower case 
letters a, b represent the numbers associated with nodes A, B, and process is described by equ. 9. 


Xf = aW C + b 
k 


(equ. 9.) 


th th 

The number within the k node of the e array is also formed in a simple manner. Since k is a 
binary number of g bits, scale it (g-e) many places to the right, filling the leftmost bits with 0's 
as they are vacated. After scaling has tak^n place, rever|g the order of the bits. This reversal, 
or inversion, produces the number in the k node of the e array. 



The rule^for locating tl^arrows of the k node in the e arra^are as follows: The solid line 

to the k node in the e array comes from the node in the e-1 array whose address (vertical 

position) i|jthe same as k's except bit k is set to 1. The dashed line originates at the node 

in the e-1 array whose address is the slime but k is set to 0. 

7 g-e 

By using the above rules, a tree-graph for any N=2® inputs can be drawn. The numbers associated 
with the nodes in the second array are computed from the inputs, using the arrow rules. Simi¬ 
larly, the numbers associated with the nodes in the third array come from the numbers associated 
with the nodes of the second array, and so on. The numbers associated with the nodes of the 
last array are actually the S(n)'s of equ. 7. 


Computationally, it is not necessary to create new arrays continually, rather, the computations 
can be done in place. This allows for a considerable saving of space, and a larger number of 
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input values. It should also be noted that the values of W can either be calculated on the spot, 
or found in a look-up table, which allows for a further saving of time. Furthermore, the creation 
of the second and third arrays, the first and second passes, can be reduced from corr^)lex multi¬ 
plications to simple additions and subtractions because of the nature and value of W . This is 
enhanced if the input data is real, with no imaginary values. 

The program that was written to perform the FFT takes samples using the analog to digital 
converter of the AX08. The sampling rate can be either fixed or variable. Display of the input 
time series is provided. The output is also in the form of a display. The final smoothed frequency 
spectrum can be viewed through a variable sized window, which is under the control of the 
parameter knobs on the front panel. The Spectrum that is displayed is a magnitude only spectrum, 
and modified in that the components a and jb are only squared and added; the square root is 
not taken, and the^t term is left out. For matters of scaling, the multiplication routine has 
been normalized so that when a number is multiplied by the largest number available in the 
computer, the answer is twice the initial. Single precision is used throughout, and in order to 
take care of fractional frequencies, and to provide statistical validity, the data was subjected 
to Hanning processes twice. The Hanning process weighs the coefficient with the sum of the 
weighed elements. The program is designed to handle any N=2^ inputs as long as the data is 
stored with the beginning of ^ imaginaries N away from the reals. Further restrictions are that 
N must be less than 1024 (=2 ), and there must be a trigonometric table for values of sine equal 

in length to the number of data points. 
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(a) Continuous function 

(b) Sampling function 
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Figure 1 

Graphical Representation of DFT 
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Figure 2 

Tree-graph for N=8 
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